################################################################################
## The code in this file is used to conduct G study and D study
################################################################################


library(gtheory)
library(matrixStats)
###############################################################################
######################      PGY2  G study   ###################################
###############################################################################
###For all levels 
PGY2 <- data.frame(matrix(ncol = 8, nrow = 12))
for (n in 1:8)
  {
  ###setwd("C:/XXX")
   dat <- read.csv(paste("PGY2_", n,".csv",sep=""))%>%
     select(-c("rater","PGY","item","X"))
   formula.dat <- "score ~ (1 | subjectID) + (1 | procID)"
   gdat <- gstudy(data = dat, formula = formula.dat, colname.strata = "subset",
                  colname.objects = "subjectID")
   ###setwd("C:/XXX")
   a <- gdat$within$`1`$components %>%
     mutate(rater="faculty") %>%
     mutate(item="item1")
   
   b <- gdat$within$`2`$components%>%
     mutate(rater="faculty") %>%
     mutate(item="item2")
   
   c <- gdat$within$`3`$components%>%
     mutate(rater="trainee") %>%
     mutate(item="item1")
   
   d <- gdat$within$`4`$components%>%
     mutate(rater="trainee") %>%
     mutate(item="item2")
   PGY2[,n] <- rbind(a,b,c,d)[,2]
}

PGY2_MSD <- cbind(rowMeans(PGY2),apply(PGY2, 1, sd),rbind(a,b,c,d)[,c(1,5,6)])
names(PGY2_MSD)[1] <- "Variance"
names(PGY2_MSD)[2] <- "SD"

write.csv(PGY2_MSD, "PGY2.csv")

###Composite 
PGY2_c <- data.frame(matrix(ncol = 3, nrow = 8))
names(PGY2_c) <- c("subject","procedure","residual")
for (n in 1:8)
{
  ###setwd("C:/XXX")
  dat <- read.csv(paste("PGY2_", n,".csv",sep=""))%>%
    select(-c("rater","PGY","item","X"))
  formula.dat <- "score ~ (1 | subjectID) + (1 | procID)"
  gdat <- gstudy(data = dat, formula = formula.dat, colname.strata = "subset",
                 colname.objects = "subjectID")
  ###setwd("C:/XXX")
  subject <- mean(gdat$within$`1`$components[1,2] +
                  gdat$within$`2`$components[1,2] +
                  gdat$within$`3`$components[1,2] +
                  gdat$within$`4`$components[1,2] +
                  2*gdat$between$var.obs[2,1] + 
                  2*gdat$between$var.obs[3,1] + 
                  2*gdat$between$var.obs[4,1] +
                  2*gdat$between$var.obs[3,2] +
                  2*gdat$between$var.obs[4,2] +
                  2*gdat$between$var.obs[4,3])
  
  proc <- mean(gdat$within$`1`$components[2,2] +
               gdat$within$`2`$components[2,2] +
               gdat$within$`3`$components[2,2] +
               gdat$within$`4`$components[2,2] +
                 2*gdat$between$var.obs[2,1] + 
                 2*gdat$between$var.obs[3,1] + 
                 2*gdat$between$var.obs[4,1] +
                 2*gdat$between$var.obs[3,2] +
                 2*gdat$between$var.obs[4,2] +
                 2*gdat$between$var.obs[4,3])
  
  resi <- mean(gdat$within$`1`$components[3,2] +
               gdat$within$`2`$components[3,2] +
               gdat$within$`3`$components[3,2] +
               gdat$within$`4`$components[3,2] +
                 2*gdat$between$var.obs[2,1] + 
                 2*gdat$between$var.obs[3,1] + 
                 2*gdat$between$var.obs[4,1] +
                 2*gdat$between$var.obs[3,2] +
                 2*gdat$between$var.obs[4,2] +
                 2*gdat$between$var.obs[4,3])
  
  
  PGY2_c[n,1] <- subject
  PGY2_c[n,2] <- proc
  PGY2_c[n,3] <- resi
}

write.csv(PGY2_c, "PGY2_c.csv")

###############################################################################
######################          PGY3        ###################################
###############################################################################
PGY3 <- data.frame(matrix(ncol = 18, nrow = 12))
for (n in 1:18)
{
  ###setwd("C:/XXX")
  dat <- read.csv(paste("PGY3_", n,".csv",sep=""))%>%
    select(-c("rater","PGY","item","X"))
  formula.dat <- "score ~ (1 | subjectID) + (1 | procID)"
  gdat <- gstudy(data = dat, formula = formula.dat, colname.strata = "subset",
                 colname.objects = "subjectID")
  ###setwd("C:/XXX")
  a <- gdat$within$`1`$components %>%
    mutate(rater="faculty") %>%
    mutate(item="item1")
  
  b <- gdat$within$`2`$components%>%
    mutate(rater="faculty") %>%
    mutate(item="item2")
  
  c <- gdat$within$`3`$components%>%
    mutate(rater="trainee") %>%
    mutate(item="item1")
  
  d <- gdat$within$`4`$components%>%
    mutate(rater="trainee") %>%
    mutate(item="item2")
  PGY3[,n] <- rbind(a,b,c,d)[,2]
}

PGY3_MSD <- cbind(rowMeans(PGY3),apply(PGY3, 1, sd),rbind(a,b,c,d)[,c(1,5,6)])
names(PGY3_MSD)[1] <- "Variance"
names(PGY3_MSD)[2] <- "SD"

write.csv(PGY3_MSD, "PGY3.csv")

###Composite 
PGY3_c <- data.frame(matrix(ncol = 3, nrow = 18))
names(PGY3_c) <- c("subject","procedure","residual")
for (n in 1:18)
{
  ###setwd("C:/XXX")
  dat <- read.csv(paste("PGY3_", n,".csv",sep=""))%>%
    select(-c("rater","PGY","item","X"))
  formula.dat <- "score ~ (1 | subjectID) + (1 | procID)"
  gdat <- gstudy(data = dat, formula = formula.dat, colname.strata = "subset",
                 colname.objects = "subjectID")
  ###setwd("C:/XXX")
  subject <- mean(gdat$within$`1`$components[1,2] +
                    gdat$within$`2`$components[1,2] +
                    gdat$within$`3`$components[1,2] +
                    gdat$within$`4`$components[1,2] +
                    gdat$between$var.obs[2,1] + 
                    gdat$between$var.obs[3,1] + 
                    gdat$between$var.obs[4,1] +
                    gdat$between$var.obs[3,2] +
                    gdat$between$var.obs[4,2] +
                    gdat$between$var.obs[4,3])
  
  proc <- mean(gdat$within$`1`$components[2,2] +
                 gdat$within$`2`$components[2,2] +
                 gdat$within$`3`$components[2,2] +
                 gdat$within$`4`$components[2,2] +
                 2*gdat$between$var.obs[2,1] + 
                 2*gdat$between$var.obs[3,1] + 
                 2*gdat$between$var.obs[4,1] +
                 2*gdat$between$var.obs[3,2] +
                 2*gdat$between$var.obs[4,2] +
                 2*gdat$between$var.obs[4,3])
  
  resi <- mean(gdat$within$`1`$components[3,2] +
                 gdat$within$`2`$components[3,2] +
                 gdat$within$`3`$components[3,2] +
                 gdat$within$`4`$components[3,2] +
                 2*gdat$between$var.obs[2,1] + 
                 2*gdat$between$var.obs[3,1] + 
                 2*gdat$between$var.obs[4,1] +
                 2*gdat$between$var.obs[3,2] +
                 2*gdat$between$var.obs[4,2] +
                 2*gdat$between$var.obs[4,3])
  
  
  PGY3_c[n,1] <- subject
  PGY3_c[n,2] <- proc
  PGY3_c[n,3] <- resi
}

write.csv(PGY3_c, "PGY3_c.csv")


          
###############################################################################
######################          PGY4        ###################################
###############################################################################
### for all levels 
PGY4 <- data.frame(matrix(ncol = 18, nrow = 12))
for (n in 1:18)
{
  ###setwd("C:/XXX")
  dat <- read.csv(paste("PGY4_", n,".csv",sep=""))%>%
    select(-c("rater","PGY","item","X"))
  formula.dat <- "score ~ (1 | subjectID) + (1 | procID)"
  gdat <- gstudy(data = dat, formula = formula.dat, colname.strata = "subset",
                 colname.objects = "subjectID")
  ###setwd("C:/XXX")
  a <- gdat$within$`1`$components %>%
    mutate(rater="faculty") %>%
    mutate(item="item1")
  
  b <- gdat$within$`2`$components%>%
    mutate(rater="faculty") %>%
    mutate(item="item2")
  
  c <- gdat$within$`3`$components%>%
    mutate(rater="trainee") %>%
    mutate(item="item1")
  
  d <- gdat$within$`4`$components%>%
    mutate(rater="trainee") %>%
    mutate(item="item2")
  PGY4[,n] <- rbind(a,b,c,d)[,2]
}

PGY4_MeanSD <- cbind(rowMeans(PGY4),apply(PGY4, 1, sd),rbind(a,b,c,d)[,c(1,5,6)])
names(PGY4_MeanSD)[1] <- "Variance"
names(PGY4_MeanSD)[2] <- "SD"

write.csv(PGY4_MeanSD, "PGY4.csv") 

### Composite
PGY4_c <- data.frame(matrix(ncol = 3, nrow = 18))
names(PGY4_c) <- c("subject","procedure","residual")
for (n in 1:18)
{
  ###setwd("C:/XXX")
  dat <- read.csv(paste("PGY4_", n,".csv",sep=""))%>%
    select(-c("rater","PGY","item","X"))
  formula.dat <- "score ~ (1 | subjectID) + (1 | procID)"
  gdat <- gstudy(data = dat, formula = formula.dat, colname.strata = "subset",
                 colname.objects = "subjectID")
  ###setwd("C:/XXX")
  subject <- mean(gdat$within$`1`$components[1,2] +
                    gdat$within$`2`$components[1,2] +
                    gdat$within$`3`$components[1,2] +
                    gdat$within$`4`$components[1,2] +
                    gdat$between$var.obs[2,1] + 
                    gdat$between$var.obs[3,1] + 
                    gdat$between$var.obs[4,1] +
                    gdat$between$var.obs[3,2] +
                    gdat$between$var.obs[4,2] +
                    gdat$between$var.obs[4,3])
  
  proc <- mean(gdat$within$`1`$components[2,2] +
                 gdat$within$`2`$components[2,2] +
                 gdat$within$`3`$components[2,2] +
                 gdat$within$`4`$components[2,2] +
                 2*gdat$between$var.obs[2,1] + 
                 2*gdat$between$var.obs[3,1] + 
                 2*gdat$between$var.obs[4,1] +
                 2*gdat$between$var.obs[3,2] +
                 2*gdat$between$var.obs[4,2] +
                 2*gdat$between$var.obs[4,3])
  
  resi <- mean(gdat$within$`1`$components[3,2] +
                 gdat$within$`2`$components[3,2] +
                 gdat$within$`3`$components[3,2] +
                 gdat$within$`4`$components[3,2] +
                 2*gdat$between$var.obs[2,1] + 
                 2*gdat$between$var.obs[3,1] + 
                 2*gdat$between$var.obs[4,1] +
                 2*gdat$between$var.obs[3,2] +
                 2*gdat$between$var.obs[4,2] +
                 2*gdat$between$var.obs[4,3])
  
  
  PGY4_c[n,1] <- subject
  PGY4_c[n,2] <- proc
  PGY4_c[n,3] <- resi
}

write.csv(PGY4_c, "PGY4_c.csv")

###############################################################################
######################          PGY5        ###################################
###############################################################################
PGY5 <- data.frame(matrix(ncol = 32, nrow = 12))
for (n in 1:32)
{
  ###setwd("C:/XXX")
  dat <- read.csv(paste("PGY5_", n,".csv",sep=""))%>%
    select(-c("rater","PGY","item","X"))
  formula.dat <- "score ~ (1 | subjectID) + (1 | procID)"
  gdat <- gstudy(data = dat, formula = formula.dat, colname.strata = "subset",
                 colname.objects = "subjectID")
  ###setwd("C:/XXX")
  a <- gdat$within$`1`$components %>%
    mutate(rater="faculty") %>%
    mutate(item="item1")
  
  b <- gdat$within$`2`$components%>%
    mutate(rater="faculty") %>%
    mutate(item="item2")
  
  c <- gdat$within$`3`$components%>%
    mutate(rater="trainee") %>%
    mutate(item="item1")
  
  d <- gdat$within$`4`$components%>%
    mutate(rater="trainee") %>%
    mutate(item="item2")
  PGY5[,n] <- rbind(a,b,c,d)[,2]
}

PGY5_MeanSD <- cbind(rowMeans(PGY5),apply(PGY5, 1, sd),rbind(a,b,c,d)[,c(1,5,6)])
names(PGY5_MeanSD)[1] <- "Variance"
names(PGY5_MeanSD)[2] <- "SD"

write.csv(PGY5_MeanSD, "PGY5.csv")    

### Composite
PGY5_c <- data.frame(matrix(ncol = 3, nrow = 32))
names(PGY5_c) <- c("subject","procedure","residual")
for (n in 1:32)
{
  ###setwd("C:/XXX")
  dat <- read.csv(paste("PGY5_", n,".csv",sep=""))%>%
    select(-c("rater","PGY","item","X"))
  formula.dat <- "score ~ (1 | subjectID) + (1 | procID)"
  gdat <- gstudy(data = dat, formula = formula.dat, colname.strata = "subset",
                 colname.objects = "subjectID")
  ###setwd("C:/XXX")
  subject <- mean(gdat$within$`1`$components[1,2] +
                    gdat$within$`2`$components[1,2] +
                    gdat$within$`3`$components[1,2] +
                    gdat$within$`4`$components[1,2] +
                    gdat$between$var.obs[2,1] + 
                    gdat$between$var.obs[3,1] + 
                    gdat$between$var.obs[4,1] +
                    gdat$between$var.obs[3,2] +
                    gdat$between$var.obs[4,2] +
                    gdat$between$var.obs[4,3])
  
  proc <- mean(gdat$within$`1`$components[2,2] +
                 gdat$within$`2`$components[2,2] +
                 gdat$within$`3`$components[2,2] +
                 gdat$within$`4`$components[2,2] +
                 2*gdat$between$var.obs[2,1] + 
                 2*gdat$between$var.obs[3,1] + 
                 2*gdat$between$var.obs[4,1] +
                 2*gdat$between$var.obs[3,2] +
                 2*gdat$between$var.obs[4,2] +
                 2*gdat$between$var.obs[4,3])
  
  resi <- mean(gdat$within$`1`$components[3,2] +
                 gdat$within$`2`$components[3,2] +
                 gdat$within$`3`$components[3,2] +
                 gdat$within$`4`$components[3,2] +
                 2*gdat$between$var.obs[2,1] + 
                 2*gdat$between$var.obs[3,1] + 
                 2*gdat$between$var.obs[4,1] +
                 2*gdat$between$var.obs[3,2] +
                 2*gdat$between$var.obs[4,2] +
                 2*gdat$between$var.obs[4,3])
  
  
  PGY5_c[n,1] <- subject
  PGY5_c[n,2] <- proc
  PGY5_c[n,3] <- resi
}

write.csv(PGY5_c, "PGY5_c.csv")


###############################################################################
######################      PGY2  D study   ###################################
###############################################################################
PGY2_D <- data.frame(matrix(ncol = 4, nrow = 8))
for (n in 1:8)
{
  ###setwd("C:/XXX")
  
  dat <- read.csv(paste("PGY2_", n,".csv",sep=""))%>%
    select(-c("rater","PGY","item","X"))
  formula.dat <- "score ~ (1 | subjectID) + (1 | procID)"
  gstudy.out <- gstudy(data = dat, formula = formula.dat,
                       colname.strata = "subset", colname.objects = "subjectID")
  ddat <- dstudy(gstudy.out, colname.objects = "subjectID", data = dat, colname.scores = "score",
         colname.strata = "subset")
  
  ###setwd("C:/XXX")
  var.universe <- ddat$composite$var.universe
  var.error.abs <- ddat$composite$var.error.abs 
  var.error.rel <- ddat$composite$var.error.rel
  gen <- ddat$composite$generalizability
  dep <- ddat$composite$dependability
  PGY2_D[n,1] <- var.universe[1]
  PGY2_D[n,2] <- var.error.abs[1]
  PGY2_D[n,3] <- var.error.rel[1]
  PGY2_D[n,4] <- gen[1]
  PGY2_D[n,5] <- dep[1]
}

names(PGY2_D)<- c("var.universe","var.error.abs", "var.error.rel","generalizability","dependability")
write.csv(PGY2_D, "PGY2_D.csv")

###############################################################################
######################      PGY3  D study   ###################################
###############################################################################
PGY3_D <- data.frame(matrix(ncol = 4, nrow = 18))
for (n in 1:18)
{
  ###setwd("C:/XXX")
  
  dat <- read.csv(paste("PGY3_", n,".csv",sep=""))%>%
    select(-c("rater","PGY","item","X"))
  formula.dat <- "score ~ (1 | subjectID) + (1 | procID)"
  gstudy.out <- gstudy(data = dat, formula = formula.dat,
                       colname.strata = "subset", colname.objects = "subjectID")
  ddat <- dstudy(gstudy.out, colname.objects = "subjectID", data = dat, colname.scores = "score",
                 colname.strata = "subset")
  
  ###setwd("C:/XXX")
  var.universe <- ddat$composite$var.universe
  var.error.abs <- ddat$composite$var.error.abs 
  var.error.rel <- ddat$composite$var.error.rel
  gen <- ddat$composite$generalizability
  dep <- ddat$composite$dependability
  PGY3_D[n,1] <- var.universe[1]
  PGY3_D[n,2] <- var.error.abs[1]
  PGY3_D[n,3] <- var.error.rel[1]
  PGY3_D[n,4] <- gen[1]
  PGY3_D[n,5] <- dep[1]
}

names(PGY3_D)<- c("var.universe","var.error.abs", "var.error.rel","generalizability","dependability")
write.csv(PGY3_D, "PGY3_D.csv")

###############################################################################
######################      PGY4  D study   ###################################
###############################################################################
PGY4_D <- data.frame(matrix(ncol = 4, nrow = 18))
for (n in 1:18)
{
  ###setwd("C:/XXX")
  
  dat <- read.csv(paste("PGY4_", n,".csv",sep=""))%>%
    select(-c("rater","PGY","item","X"))
  formula.dat <- "score ~ (1 | subjectID) + (1 | procID)"
  gstudy.out <- gstudy(data = dat, formula = formula.dat,
                       colname.strata = "subset", colname.objects = "subjectID")
  ddat <- dstudy(gstudy.out, colname.objects = "subjectID", data = dat, colname.scores = "score",
                 colname.strata = "subset")
  
  ###setwd("C:/XXX")
  var.universe <- ddat$composite$var.universe
  var.error.abs <- ddat$composite$var.error.abs 
  var.error.rel <- ddat$composite$var.error.rel
  gen <- ddat$composite$generalizability
  dep <- ddat$composite$dependability
  PGY4_D[n,1] <- var.universe[1]
  PGY4_D[n,2] <- var.error.abs[1]
  PGY4_D[n,3] <- var.error.rel[1]
  PGY4_D[n,4] <- gen[1]
  PGY4_D[n,5] <- dep[1]
}
names(PGY4_D)<- c("var.universe","var.error.abs", "var.error.rel","generalizability","dependability")
write.csv(PGY4_D, "PGY4_D.csv")

###############################################################################
######################      PGY5 D study   ###################################
###############################################################################
PGY5_D <- data.frame(matrix(ncol = 4, nrow = 32))
for (n in 1:32)
{
  ###setwd("C:/XXX")
  
  dat <- read.csv(paste("PGY5_", n,".csv",sep=""))%>%
    select(-c("rater","PGY","item","X"))
  formula.dat <- "score ~ (1 | subjectID) + (1 | procID)"
  gstudy.out <- gstudy(data = dat, formula = formula.dat,
                       colname.strata = "subset", colname.objects = "subjectID")
  ddat <- dstudy(gstudy.out, colname.objects = "subjectID", data = dat, colname.scores = "score",
                 colname.strata = "subset")
  
  ###setwd("C:/XXX")
  general <- ddat$composite$generalizability 
  depend <- ddat$composite$dependability
  var.universe <- ddat$composite$var.universe
  sem.abs <- ddat$composite$sem.abs
  var.universe <- ddat$composite$var.universe
  var.error.abs <- ddat$composite$var.error.abs 
  var.error.rel <- ddat$composite$var.error.rel
  gen <- ddat$composite$generalizability
  dep <- ddat$composite$dependability
  PGY5_D[n,1] <- var.universe[1]
  PGY5_D[n,2] <- var.error.abs[1]
  PGY5_D[n,3] <- var.error.rel[1]
  PGY5_D[n,4] <- gen[1]
  PGY5_D[n,5] <- dep[1]
}

names(PGY5_D)<- c("var.universe","var.error.abs", "var.error.rel","generalizability","dependability")
write.csv(PGY5_D, "PGY5_D.csv")








